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differential equations with 
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Abstract 


The class of strong stochastic Runge-Kutta (SRK) methods for stochas- 
tic differential equations with a commutative noise proposed by Rofler (2010) 
is considered. Motivated by Komori and Burrage (2013), we design a class of 
explicit stochastic orthogonal Runge-Kutta Chebyshev (SROCKC2) meth- 
ods of strong order one for the approximation of the solution of Its SDEs 
with an m-dimensional commutative noise. The mean-square and asymptotic 
stability analysis of the newly proposed methods applied to a scalar linear 
test equation with a multiplicative noise is presented. Finally, some numer- 
ical experiments for stochastic models arising in applications are given that 
confirm the theoretical discussion. 


Keywords: Stochastic differential equations; Runge-Kutta methods; Stochas- 
tic mean square stability, Stiff equations; Commutative noise. 


1 Introduction 


Stochastic differential equations (SDE s) arise in a variety of applied mathe- 
matics and physics areas; see [4,22]. Such equations are the results of adding 
random fluctuations in the parameters of a system of ordinary differential 
equations (ODEs). In the literature, there are plenty of work devoted to the 
study of stochastic differential calculus; see [10,12]. However, in many cases, 
the analytical solutions of such SDEs are not known; this makes numerical 
methods as very important tools for solving SDEs, which have been consid- 
ered by many researchers. For solving different forms of SDEs, many efficient 
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numerical methods have been constructed, for example, [6, 12, 16,23]. Espe- 
cially, SRK methods are an important derivative free class of numerical meth- 
ods, which are proposed for strong approximations of SDEs; see [5, 15, 17]. 
In 2010, RoBler [19,20] introduced a new class of SRK methods in which the 
number of stages does not depend on the dimension, m, of the Wiener pro- 
cess. On the other hand, the stability of numerical methods for SDEs is an 
essential factor in order to avoid possible explosion. Generally, the implicit 
methods have a wider range of acceptable step sizes compared to explicit ones, 
which makes them suitable for the solution of stiff systems; see [7,8, 12, 24]. 
However, the implementation of implicit methods demands to solve a nonlin- 
ear system of equations per step, which may increase the computational cost 
significantly in the case of high dimensional SDEs. Therefore, explicit meth- 
ods with extended stability regions could greatly help to solve stiff SDEs. In 
this regard, Abdulle et al. proposed two classes of explicit stochastic orthog- 
onal Runge-Kutta Chebyshev (SROCK) methods of strong order 1/2 for 
SDEs in Stratonovich [1] and Ité [2] sense. Recently, Komori et al. in [14], 
constructed explicit strong SRK methods (SROCK D1 and SROCK D2) of 
order one for SDEs in It6 and Stratonovich sense. These methods are de- 
signed for approximation of SDEs with noncommutative noises. As we know, 
double stochastic integrals must be used in methods of strong order one for 
the approximation of the solution of a SDE with a noncommutative noise; 
see [12,13]. Such an expensive computational cost can be reduced when the 
SDE has a single Wiener process or has a commutative noise. It should be 
noted, such SDEs can be used for modeling a variety of nonlinear physical 
systems; see [12]. In this article, the SRK methods of [20] for the strong 
approximation of SDEs with a commutative noise is briefly reviewed. Then, 
motivated by [14], a new class of SROCK methods of strong order one for 
the approximation of the solution of SDEs with an m-dimensional commu- 
tative noise is proposed. The mean-square stability (MS-stability) functions 
for the new proposed methods, applied to a scalar linear test equations with 
a multiplicative noise, are determined, and in some special cases, the regions 
of MS-stability and T-stability are compared between our proposed methods 
and the scalar test equation. 


Consider the autonomous SDE of the It6 type in the form of 
dX(t) =a(X(t))dt + 0(X(t))aW, X(to)=Xo, te lto,T] (1) 


Here, X(t) € R4, a: [to, T] x R¢ > R%, is a drift vector, b = (bj, b2,...,bm): 
[to, T] x R? > R?@*" is a diffusion matrix, and W(t) is an m-dimensional 
standard Wiener process defined on probability space (Q, 7, P) with the fil- 
tration {F¢}+e[to,r], Which satisfies the usual conditions. It is assumed, Xo is 
F;,, -measurable and independent of the Wiener process with E|Xo|? < oo, 
where |.| denotes the Euclidean norm. Suppose that the conditions of the ex- 
istence and uniqueness theorem [12] are fulfilled for SDE (1). For simplicity 
in this paper, the equidistant discretization py, = {to,t1,...,t~ = T} of the 
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time interval [to,T] with t, = to + nh, where h = ote for n = 0,1,...,N, 
has been considered. Also, for an equidistant grid point t,, € pp, the notation 
Yn denotes a discrete approximation to the solution X(t,,) of (1). 

In the following, the SDE (1) is considered with a commutative noise, in 
the other words, the diffusion matrix b satisfies the commutativity condition 


d 
. per 
4=1 


fot alle = 1am, k= 1,040, Sad (2) -€ lip, T | eR 


ov 


ayes = gg0br 
aor (2) =) (t,2), (2) 
i=l 


et 


Definition 1 (see [12,20]). A discrete time approximation y” is said to be 
convergent with the strong order a to the solution X of SDE (1) at time tn, 
if there exist constants C > 0 and do > 0 such that, for each h € (0,60), 


cE 
(E1X (tn) - yh?) " < Che (3) 
The constants C' and 69 are independent of h. 


The outline of the paper is as follows: In Section 2, we will introduce 
second order Chebyshev methods proposed by Abdulle and Medovikov in [3] 
for ODEs. In Section 3, we introduce the class of SRK methods proposed 
by Rofler in [20] for the strong approximation of SDEs with a commutative 
noise. Then, we introduce a new class of methods of strong order one for the 
approximation of the solution of It6 SDEs with an m-dimensional commuta- 
tive noise. In Section 4, a brief overview of the stochastic stability concepts of 
SDEs is presented. We also analyze the MS-stability and T-stability proper- 
ties of the newly proposed method on a linear scalar SDE. In Section 5, some 
numerical results for stochastic models arising in applications are provided 
to confirm theoretical discussions. 


2 Chebyshev methods of order two for ODEs 


Consider the autonomous d-dimensional ODE system X‘(t) = a(X(t)) with 
the initial condition X(to) = yo. The s-stage explicit deterministic Runge— 
Kutta (RK) method for solving this system corresponding to step-size h is 
as follows: 


Yntt =Yn thd. aa(H), (4) 


where 
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Applying (4) to the scalar test equation yields 
X(t) =AX(t), t>0, y(0) = yo, (5) 


where R(A) < 0 and yo ¥ 0, we have ynii = Rs(Ah)yYn. Here, R,(z) is 
called a stability function. Consequently, the stability region of (4) can be 
expressed as 

Rs = {z €C: |Rs(z)| < 1}. (6) 


Riha [18] showed that the polynomial 


2 


s 
R(z)=1l+z2+ 2 +2 Csj ER, 
jJ= 


exists such that |R,(z)| < 1 for z € [—1,,0] with J, as large as possible. It 
is desirable in practice to replace |R,(z)| < 1 by |R,(z)| <7 < 1 (damping). 
Abdulle et al. in [3] constructed approximations to these optimal stability 
polynomials. In this approach, at first for a given 7 (damping factor) and 
s (stage value), the values a, and 1, and also the polynomial w(x) of the 
degree two with complex zeros a, + 3, which satisfies w(a,) = 1, are cal- 
culated. Then, the orthogonal polynomials {Q;(z) = with respect to the 


weight function le on [—1,1] are constructed in which Q;(a;) = 1 for 


j =0,1,...,s—2. Finally, the approximation of optimal stability polynomi- 
als is constructed as 


Rs(x) = w(x)Qs—2(2), (7) 
where w(x) := W(as+2/d,) with ds = te > 0 and Q;(z) := Q;(as+2/ds) 
for 7 = 0,1,...,s— 2. For more details, you can refer to [3]. The polyno- 
mials Qo(x), Q1(2),...,Qs—2(x) satisfy in the below three-term recurrence 
relations as form 

Qo(z)=1, Qi(z)=14+ mz, (8) 
Q5 (2) = (Hjz + Ky + 1)Q5-1(2) — Kj Q5-2(z), 
for 7 =1,...,s—2. The values of ; and «; can be determined by inserting 


the different suitable values 71, r2 € R\ {0} in (8) and solving the nonsingular 
linear system 


(uri + &j + 1)Q5-1(ri) — &jQ5-2(ri) = Q3(ri), 7= 1,2. 


Abdulle et al. in [3], corresponding to stability polynomial (7), con- 
structed the second order explicit RK method {y,},>0 for the approximation 
of the solution of the deterministic system X'(t) = a(X(t)). The mentioned 
method is implemented in a code called ROCK2 as below: 


0 
rel ) = Un; 
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HO = = Yn + hya(H}”), 


HS, = bys (A SY + (ip + 1H” — eH, 2<j<s-2, 
H®) = HO, + ho,a(H), 

Ts 0 
ner = HY? + hOsa(HY) — h4(1 — ge )(a(H) — a(H,)), (9) 


s 


where the parameters 0, and T, satisfy in the relation w(2) = 1+20,7+7,27. 


3 SRK methods for SDEs with a commutative noise 


In this section, we will review the SRK method introduced by Roéfler in [20] 
for the strong approximation of Ito SDEs with m-dimensional commutative 
noise. In this class, the SRK method is given by yo = Xo and 


Yn4+1 = Yn + S- aja(ty + OO An, HO) hn 
i=l 


+50 578 Tag Bale Gite lin HO); 0) 


i=1 k=1 


for n = 0,1,2,..., with stage-values 


HO = yt So AM alty + hn, HO Yan + >> BOB (ty + Oh, HOY 0, 


= j=l l=1 
H® = yn, + SAD alt, + ty HO) )iig — 57 Be BF (ty +e hy HP) Vin 
j=l j=l 
+ ‘2 > BOO (ty + oY hy, HY?) as 
j=l I=1 
fori =1,...,s and k =1,...,m. The random variables of the method are 


defined by the stochastic Ito integrals 


tn th 
Ie) n = W* (tn +h) — W* (tn) = dw®. 
ty 
Since the random variables of the method (10) are only increments of the 
Wiener process and the simulation of the multiple stochastic integrals is not 
required any more, the computational cost will be reduced significantly. 
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Table 1: Butcher tableau for the SRK methods SRIC1 (left) and SRIC2 (right) 


0 0 
0] 0 0 0} 1 0 
0|0 0 0 0 0|0 0 0 0 
0 0 
0} 0 1 Lid 1 
0/0 0 |-l 2 oe ee 

1 0 O0[1 0 0 3 5 Of 1 0 0 

0s 3 ae 


In the following, let C?*4([to, T] x R¢,R“) denote the space of all mappings 
from [to,T] x R¢ to R¢, which are p and q times continuously differentiable 
with respect to the time and the state variables, respectively. 


Remark 1. Frequently, the coefficients of the SRK method (10) are shown 
by the following Butcher tableau : 


ce) | 4 | BO 

ay | Am | Bw 

at | go" 
p27 


Theorem 1 (Order conditions). Let a,b’ € C!?([to,T] x R?,R) for j = 
1,...,m. Then, the SRK method (10) obtains the strong order 0.5 for ap- 
proximation of the solution of (1), if its coefficients satisfy in the equations 


1. afe=1, 2 BY e=1, 3. B@e=0. 


In addition, suppose a,b) € Cl?([to,T] x R4,R) for j =1,...,m. Then, the 
SRK method (10) obtains the strong order 1.0, if the equations 


4, B2”ADe=0, 5. BQ (BOBMe) =0, 
6 BY" BYe=0 7. B@"(BMe)2=0, 8. BO” BMe =1, 


are fulfilled with C = AMe for i= 1,2. 


Proof. see [20]. 


To date, suggestions have been presented for the coefficients of SRK 
method (10). As an example, R6fler in [20] purposed three stages of 
SRK methods SRIC1 and SRIC2 with orders (pp,ps) = (1.0,1.0) and 
(pp, ps) = (2.0,1.0), respectively, in which pp is the order of the method 
related to ODEs, that is, b= 0 and pg is the order of the method related to 
SDEs; see Table 1. 
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Remark 2. In the rest of the paper, we set B©) =0 and 


0--- 0 00 D . 

BM =!o.. 000) -e=]2] .o%=] 2] an 
0.» 100 - ; 
O—100 2, 

sxs 0 eee 
sxl 2 sxl1 


therefore, the equations 2. and 3. as well as the equations 5.—8. in Theorem 
1 are satisfied. 


3.1 A new class of Chebyshev methods for SDEs 


In this section, motivated by Komori et al. [14], we introduce a new class of 
s-stage explicit stochastic orthogonal Runge-Kutta Chebyshev method for 
Ito SDEs with a commutative noise. We denote this method by SROCKC2, 
which has the orders (pp, pg) = (2.0, 1.0). 


For the strong approximation of the solution of SDE (1) with a commu- 
tative noise, the s-stage (s > 2) SROCKC2 method is given by yo = Xo 
and 


yn = 1) + 2 (BO Tn + BO? Vin) O* (tn + rn, HO), (12) 


i=1 k=1 
with stage values 

He” = Yn, HY = Yn + hura(H{), 

HA, = hyjga(HS”) + (mj + HY — 6 HO, 2<j<s-2, 


HS) = HY By hd,a(H.,), 
Hy), = HO + hO.a(HQ) — hO6(1 — 7)(a( HL) - a( H%)), 


s 


HY =y,, + s AY a(t ha = 3 BD Hm) 
j=l j=l 
+d apie) Baton, 
j=l l=1 me 
where pj; and Kk; (j =1,...,5— 2) as well as @, and 7, are free parameters 


in R, which are in accordance with the second method ROCK2 in (9); see [3] 
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for more details. In addition, B®, B™, 6 and B®) would be defined as 
previously mentioned in relation (11). 


Remark 3. We note that in the deterministic case, b = 0, the method (12) 
coincides with the second order ROCK2 method, introduced in (9). There- 
fore, assuming Yn + >>;_4 aja(H!) Es H®, from second order properties of 
the ROCK2 method, we have S>?_, a; = 1, and the equation 1. in Theorem 
1 is satisfied. 


A Stability properties 


In this section, we will first investigate the MS-stability of our proposed 
methods and compare them with the MS-stability regions of the methods 
SRIC1 and SRIC2 in [20]. Then, the T-stability properties of the proposed 
schemes on a scalar linear test equation with a multiplicative noise have 
been analyzed. In this regard, we first review the MS-stability analysis of the 
methods SRIC1 and SRIC2 on the linear test problem with the multiplicative 
noise 

dX (t) = AX (t)dt + uX(t)dW;, A, WE C, (13) 


with the deterministic initial condition X (to) = yo € R\{0}. For A,y € C, 
the solution of (13), X(t) = yoexp{(A — $u*)(t — to) + u(Wi — Wi,)}, is 
stochastically asymptotically stable in the large if and only if 


limy-00 |X (t)| = 0 R(A — 5p?) < 0. (14) 


Also, the equilibrium position of SDE (13) is asymptotically MS-stable if and 
only if 
limt+oo E(|X(#)|?) = 0 & 2R(A) + |u|? < 0, (15) 


for A, € C; see [12]. It worth mentioning, MS-stability always results in 
asymptotically stability in the large due to inequality (2A — yu?) < 2R(A) + 
|u|?. 

Here, we are going to find conditions in which the SRIC1 and SRIC2 
methods, on the scalar test equation (13), have numerically stable solutions. 
For every A, € C a numerical method applied on the test equation (13) is 
numerically asymptotically stable and MS-stable, respectively, if and only if 
limp—sco |Yn| = 0 (with probability 1) and limy_.. E(|yn|?) = 0. To deter- 
mine the domain of asymptotic stability of a numerical method, we recall the 
following lemma from [9]. 


Lemma 1. For a given sequence of real-valued, nonnegative, independent, 
and identically distributed random variables {R,(h, A, L)yn}9- Consider 
the sequence of random variable {\yn|}?<o defined by yn41 = Rn(h, A, Wyn 
and |yo| 4 0 (w.p.1). If the random variables log(|Rn(h, A, 4)|) are square 
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integrable, then 


limy-+oo [Yn] = 02 B(log(|Rn(h, A, 1)|) <0. (16) 


So the domain of asymptotic stability of a numerical method can be ex- 
pressed as 


Rag pete B(log(|Rn(h, ,#)I) <0}. 


Set a, := ie and ag := l+a+ ig? ¥, where x = AX and 
y = pVvh. Applying SRIC1 and SRIC2 methods on the scalar test equation 
(13), respectively, we have 


Ta),n Tey ,m 2 . 
ner = (a; + ty + Sy on, 9 = 1,2. 1 
Yn+1 (a, Vat on 4 )y j (17) 


Here, we find the forms 


Yn+1 = Ri (x, y)Yn: J = 1, 2, 


for (17), and the MS-stability functions of the methods SRIC1 and SRIC2 
will be given, respectively, by 


E(\ynsil?) = RA(@,y) Eyl”), § = 1,2. 


Obviously, for A, 4 € C and h > 0, the methods SRIC1 and SRIC2 are MS- 
stable, if R?, (x,y) < 1 for 7 = 1,2. So, the domain of MS-stability of SRIC1 
and SRIC2 methods can be expressed, respectively, as 


Rigg = {A,HEC: Bi (a,y) <1}, 9 =1,2 


Because of better visualization of the regions Rj,, and R%,, in the x — y 
plane, the restriction 4,4 € R is considered. So, one can express the MS- 
stability function of the SRIC1 and SRIC2 methods, respectively, as 


- ° 
Ri (x,y) =apt+y?+ gy F945, 7 = 1,2. 


It should be noted that all the following figures, demonstrating regions of MS- 
stability for the numerical methods under consideration, are plotted, using 
the Mathematica software. Figure 1 gives the MS-stability regions of the 
methods SRIC1 and SRIC2. The light gray area shows the MS-stability 
region of the test equation (13) and the dark gray area shows the MS-stability 
regions of the methods. Comparison of the illustrated regions in Figure 1 
shows that the SRIC1 and SRIC2 methods have small regions of MS-stability. 
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Figure 1: MS-stability region for the test equation (13) (light grey surface) and SRK 
methods (dark grey surface) (left) SRIC1, (right) SRIC2. 


4.1 M<S-stability analysis of SROCKC2 


Consider the method (12), for simplification similar to [14], we set all com- 
ponents of A“) as zero except 


AD AD . AD A as 1<j< s—2, 


s—2,jJ s—1,j 85) 


where A() is an s xs strictly lower triangular matrix which could be obtained 
from the following equation: 


isl 
a = Yn t SA a) 
j=l 


for i = 1,...,s. Clearly, from the definition of AM, the equation 4. in 
Theorem | is satisfied. 

Now, suppose that R,(#) = w(x)Q,_2(x) is the stability polynomial of 
the s-stage deterministic method ROCKD2 (9) with parameter 1; see [3]. 
When the method (12) is applied to the test equation (13), we find the forms 
Yn41 = Rs(h, A, wp, AWn) yn, where 


1 1 
R,(h, ,u, AWn) = (wan) + AWnu + 5[AWal?n? - shu?) Qs—2(Ah). 


Applying the expectation operator and according to the properties of 
standard normal distribution in the terms of x and y, we have 
4 


B(|R.(2,y)P) = (w(2)? +9? + 2) [Q.-2(a))?. 


Therefore, the mean square stability regions of the methods SROCKC2 is 


a=3, 4-08 s=6, 1 =0.375 


to Excessive Length + 


20F 


IS5¢ 


> 10f 


-50 -40 —30 -20 -10 0 


Figure 2: Profile of the MS-stability regions of the SROCKC2 schemes for some special 
values of 7 and s. 


RMS = {(2,y) ER? : E||R(@,y)?] <1}. (18) 


For more clear results concerning the MS-stability regions of SROCKC2 
schemes, in the following, we confine our investigation to some optimal values 
of 7 and s, which have been suggested by [14]. The regions of MS-stability 
for SROCKC2 schemes with different values of 7 and s (dark gray area) are 
illustrated in Figure 2. We can see in Figure 2 that the MS-stability region of 
methods SROCKC2 is larger than the methods SRIC1 and SRIC2. In what 
follows, SROCKC2 scheme with three stages (7 = 0.4), six stages (7) = 0.375), 
and nine stages (7 = 0.3) are denoted by SROCKC2 — 3, SROCKC2 — 6, 
and SROCKC2 — 9, respectively. 


4.2 T-stability analysis of SROCKC2 


To measure the asymptotic stability of a numerical method with respect to 
the driving process, Saito and Mitsui established the definition of T-stability 
in [21]. 
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s=3, n=0.4 s=6, 7=0.375 


Figure 3: Profile of the T-stability regions of the SROCKC2 schemes for some optimal 
values of 7 and s. 


Definition 2. Assume that the test equation (13) is stochastically asymp- 
totically stable in the large. The numerical scheme equipped with a specified 
driving process said to be T-stable if limy+oo|Yn| = 0 holds for the driving 
process. 


It can be deduced that the criterion of T-stability depends both on the 
scheme and the driving process. Accordingly, similar to the approach, which 
is used by Lemma 1, the criterion of T-stability of a numerical method with 
respect to normal random variables can express as 


LQ 12 
log T(h, A, 1) = =z | (log(|n(h, X, w)len 2? dz) <0. (19) 


It seems it is not so simple to find the closed form of (19). Thus, we use 
the MATLAB software to approximate the integral in (19). The integral 
interval (—0o0o,-++00) is approximated by [—50,50] because of the magnitude 
of the integrand in (19) becomes sufficiently small when |z| > 50. In Figure 
3, the regions of T-stability for SROCKC2 schemes with different values of 
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Figure 4: Profile of the T-stability regions of the methods SRIC1 (left) and SRIC2 (right). 


71 and s (dark gray area) have been illustrated. Similarly, the regions of 
T-stability for the methods SRIC1, SRIC2, Milstein and Platen (light gray 
area) are depicted in Figures 4 and 5. Comparison of the illustrated regions 
confirms that the proposed methods have reasonably larger regions of both 
MS-stability and T-stability. 


5 Numerical experiments 


In this section, we apply some optimal cases of SROCKC2 schemes (12) on 
some examples of It6 stochastic differential equations that arising in applica- 
tions. These examples illustrate the efficiency of the new proposed methods. 
Denote yo and X(t) as the numerical solutions and the exact solution 
at step point ty in ith simulation, respectively. The root mean square er- 
ror (RMSE) is used by simulating 1000 independent trajectories for a given 
step-size h as bellow: 


1 


Tl 1000 i) 2 2 
MsE = | —— x yO!) 
RMS (i > (tw) — YN 


i=l 


We also investigate computational costs for a given h to measure accuracy 
of proposed methods. In simulation results, the number of evaluations of 
the drift function a, of the diffusion function b, and the number of random 
numbers to be generated in every step is taken as measure of computational 
costs (S,). In the following, we will compare the methods SROCKC2-3, 
SROCKC2-6, and SROCKC2-9 with some popular methods, including Euler- 
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5 5 
4 4 
3 3 
~ > 
2 2 
1 1 
0 0 
“5 4 3 2 -1 0 5 -4 -3 -2 -1 0 
x x 


Figure 5: Profile of the T-stability regions of the methods Milstein (left) and Platen 
(right). 


Table 2: A least squares fit for the parameters C' and a, for problem 1 with A= 1, =1 
Xty =1 


SROCKC2-3 = SROCKC2-6 SROCKC2-9  Platen Milstein =SRIC1 


a 1.0514 1.0597 1.0615 1.0207 1.0109 1.0109 
log(C) 1.8067 1.7033 1.6803 2.0750 2.2317 2.2317 


Maruyama [12], Milstein [16], SRIC1 and SRIC2 [20], and also the method 
proposed by Platen [12]. 


Problem 1: The first problem is the scalar test equation (13) that is 
considered on I = [0,1] with the initial condition X,, = 1. For this setting, 
the simulation results are presented in Figure 6. In this figure, RSME versus 
step size of the methods are illustrated at T = 1. It is clearly seen in Figure 
6 that the slopes of curves appear to match well and (3) for a = 1 is valid. 
Further, with assumption RSME =~ Ch® for some constants C’ and a, we 
can write 


log(RSME) = log(C) + alog(h). (20) 


In Table 2, a least squares fit for the parameters C' and a are calculated, 
respectively. These results confirm that the proposed methods have a con- 
vergence with order 1 in mean square sense. Also, in Table 3, the RMSE for 
time steps h = 2~°,...,27° are computed. 

As a further aspect of providing comparisons of the computational efficiency 
of the proposed methods, ratio CPU times of the methods SROCKC2, Mil- 
stein, Platen, and SRIC2 to the running time of the method SRIC1 are cal- 
culated in Table 4. It can be deduced from Table 4 that the CPU time of the 


10 T 
»» Order 1 
SROCKC2-3 
—*-— - SROCKC2-6 
—— SROCKC2-9 
~ — ~ Platen 
—*— Mil 
— © — RungeC1 
f 10° : 
rt 
o 
D 
oO 
g 
oO 
o 
E 
& 10° 4 
-3 
10 : 
10° 10° 10" 
At 
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Table 3: Means of absolute errors, RMSE, for problem 1 
h SROCKC2-3 SROCKC2-3 SROCKC2-3 Platen Milstein SRIC2 
g-8 0.150 0.131 0.127 0.223 0.271 0.271 
2-6 0.082 0.072 0.070 0.121 0.145 0.145 
are 0.037 0.032 0.031 0.056 0.069 0.068 
2-8 0.017 0.015 0.014 0.027 0.035 0.034 
2-9 0.008 0.007 0.007 0.013 0.017 0.016 


method SROCKC2-3 and SRIC2 are approximately similar to the method 
SRIC1. 


In following, because we do not know exact solutions to problems, refer- 
ence solutions have been simulated with the Milstein scheme [16] based on 
the step size h = 2718. 

Problem 2: The second problem is a chemical reaction model [12]. The 
mathematical model is 


dx, (t) = (cz(t) — cgay (t)(a1(t) — 1) + 2cswa(t)) ) dt 
+21 (t) (cia) + a2dW?) 

dia(t) = (Sule — 1) — e3x9(t) — cawa(t))) dt 
4+-29(t) (s.aw? ra BadW?) 


The system parameters take the values cy = cg = 200, c3 = 100, and 
c4 = 100 while the stochastic coefficients are a, = ag = 5 and $, = By = 0.5. 
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Table 4: Ratio CPU times of the proposed methods to the running time of the method 


SRIC1, for problem 1 with A=1,u=1 Xig =lat T=1 

h SROCKC2-3. SROCKC2-6 SROCKC2-9  Platen Milstein SRIC2 
2° 0.842 1.736 2.157 0.115 0.105 0.473 
a> 1.346 1.846 3.891 0.269 0.262 1.000 
2°? 0.957 2.015 3.723 0.414 0.106 1.063 
ore 0.970 2.000 3.303 0.284 0.176 0.961 
2-3 0.984 2.035 3.533 0.329 0.203 1.011 
Pa 1.032 2.060 3.486 0.319 0.151 1.001 
or8 0.995 1.949 3.212 0.386 0.239 1.022 


Table 5: A least squares fit for the parameters C and a, for problem 2 


SROCKC2-3. SROCKC2-6 = SROCKC2-9 Platen Milstein SRIC1 
a 0.7247 0.7237 0.6993 1.5669 1.4207 0.8303 
log(C) 8.5043 8.3181 7.9808 19.5330 17.8396 10.1329 


Table 6: Means of absolute errors, RMSE, for problem 2 


h SROCKC2-3 SROCKC2-3 SROCKC2-3 SRIC2 Milstein Euler Platen 
0.01 x 27° unst. 1.17 1.21 unst. unst. unst. unst. 
0.01 x 271° 0.86 0.91 0.92 unst. unst. unst.  unst. 
0.01 x 2711 0.71 0.60 0.58 0.98 1.72 unst. 1.58 
0.01 x 2712 0.41 0.34 0.33 0.55 0.51 0.51 0.41 
0.01 x 2718 0.26 0.22 0.22 0.31 0.24 0.24 0.18 


Table 7: Ratio CPU times of the proposed methods to the running time of the method 


SRIC1, for problem 2 


h SROCKC2-3 


SROCKC2-6 


SROCKC2-9 


Platen 


Milstein SRIC2 


0.01 x 276 1.066 
0.01 x 2-7 0.935 
0.01 x 278 1.034 
0.01 x 279 1.008 
0.01 x 2719 = 0.991 
0.01 x 2714 1.006 
0.01 x 2712 1.004 
0.01 x 2743 1.007 


1.333 
1.225 
1.310 
1.288 
1.287 
1.319 
1.304 
1.314 


1.600 
1.581 
1.672 
1.661 
1.665 
1.696 
1.685 
1.684 


0.733 0.667 0.933 
0.709 0.645 0.967 
0.706 0.689 1.034 
0.711 0.678 1.016 
0.707 0.671 1.042 
0.724 0.698 1.002 
0.719 0.691 1.000 
0.730 0.699 1.009 
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Figure 7: RMSE versus step size (left) and RMSE versus computational cost (right) with 
double logarithmic (ld) measures for test problem 2. 


Table 8: A least squares fit for the parameters C and a, for problem 3 


SROCKC2-3 SROCKC2-6 SROCKC2-9  Platen Milstein SRIC1 


a 0.5454 0.5715 0.5244 1.4623 1.3736 1.4243 
log(C) 0.6772 0.8068 0.5372 6.7938 6.2749 6.5723 


The initial conditions are 71(0) = 1000 and x2(0) = 100, and the integration 
is performed on the interval [0,0.01]. The results are indicated in Figure 7 
and Tables 5-7. Table 5 shows that although all mentioned methods have 
same slope values, a, but the values of C, for the methods SROCKC2-3, 
SROCKC2-6, and SROCKC2-9 are much less that other methods. Also, from 
Table 6, it can be deduced that the methods SROCKC2-6 and SROCKC2-9 
are MS-stable for step sizes h < 7.81x 107° and h < 1.561074. Clearly, the 
simulation results suggest proposed methods as the most efficient methods 
with respect to the root mean-square errors and computational costs. Also, 
the evolution of the x;(t) and x(t) is given in Figure 8. 


Problem 3: The third problem is Marine bacteriophage infection 
model, which is a stiff SDEs; see [11]. The mathematical model is 


ds(t) = (as(t)(1 — (i(t) + s(t))) — s(t)p(t) ) at +04 (s(t) = s*)dW}, 
di(t) = (s()p(¢) = ti(t)) dt 409 (i(t) - i*)aw?, (21) 


dp(t) = (—s(¢)p(6) — p(t) + bei(t)) dt — 03 (p(t) — p* aw, 
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40 individual paths 200 40 individual paths 
mean of 1000 paths mean of 1000 paths 


0 1 2 9 | 2 
t x 1074 t x10" 
Figure 8: Problem 2: The functions 21(t) and a2(t) averaged over 1000 discretized 
Brownian paths and along 40 individual paths with the SROCKC2-6 method with h = 
0.01 x 2-8. 


Table 9: Means of absolute errors, RMSE, for problem 3 


h SROCKC2-3 SROCKC2-3 SROCKC2-3 SRIC1 Milstein Euler Platen 
5x 2-5 unst. 0.353 0.351 unst. unst. unst.  unst. 
5 x 276 0.354 0.350 0.344 unst. unst. unst.  unst. 
5x2-7 0.312 0.311 0.309 unst. unst. unst. unst. 
5x 2-8 0.238 0.237 0.235 unst. unst. unst.  unst. 
5 x 27-9 0.162 0.161 0.159 unst. unst. unst. unst. 


5x 2-% 0.108 0.107 0.105 0.365 0.355 0.362 0.372 
ee ie 0.074 0.072 0.073 0.136 0.137 0.133 0.135 
Seo 0.055 0.054 0.055 0.039 0.039 0.038 0.039 


with the initial condition (s0,%9,po) = (0.3,0.2,5). Here s represents the 
susceptible bacteria, i is the infected bacteria, and p is the phage (viruses). 
The drift coefficients have the values a = 8.65, 0 = 24.628, m = 14.925, and 
b = 60 while the noise coefficients have the values 0, = 02 = o3 = 0.4 and 


as*(1—s*) ,  af(1—s*) 


Ss = = , = 
é+as* a é+as* 


In Figure 9, we report the errors versus step size of the methods and errors 
versus computational effort with double logarithmic (1d) measures on the in- 
terval [0,5]. Also, Figure 9 illustrates the reduction in computational effort 
of methods SROCKC2-3, SROCKC2-6, and SROCKC2-9 in the same level of 
accuracy. Clearly, the simulation results in Tables 8-10 suggest the proposed 
methods as the most efficient methods with respect to root mean-square er- 
rors and computational costs. The evolution of the three interacting species 
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Figure 9: RMSE versus step size (left) and RMSEs versus computational cost (right) 
with double logarithmic (ld) measures for test problem 3. 


Table 10: Ratio CPU times of the proposed methods to the running time of the method 
SRIC1, for problem 3 


h SROCKC2-3. SROCKC2-6 SROCKC2-9 Platen Milstein SRIC2 
5x2-> 1.200 1.600 3.000 0.400 0.401 1.200 
5x 278 0.909 1.545 2.454 0.363 0.181 1.000 
5x 27-7 1.001 1.476 2.428 0.381 0.285 1.095 
5x 2-8 1.024 1.561 2.243 0.390 0.268 1.048 
5x 279 0.964 1.511 2.178 0.392 0.272 1.023 
5 x 2710 1.189 1.561 2.298 0.445 0.274 1.158 
5x27 0.972 1.413 2.011 0.356 0.263 1.087 
5x4 0.952 1.384 2.158 0.371 0.267 1.069 


is given in Figure 10. 


6 Conclusions 


In this paper, the class of strong stochastic Runge-Kutta methods for stochas- 
tic differential equations with commutative noise due to Rofler (2010) is 
considered. For stiff SDEs, a family of explicit stochastic orthogonal Runge— 
Kutta Chebyshev methods of strong order one for the approximation of the 
solution of It6 SDEs with m-dimensional commutative noise are designed. 
The mean-square and asymptotic stability of newly proposed methods ap- 
plied to a scalar linear test equation with multiplicative noise have been an- 
alyzed. Some numerical results for stochastic models arising in applications 
are provided to confirm theoretical discussions. In future works, we will con- 
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Figure 10: Problem 3: The functions s(t), i(t), and p(t) averaged over 1000 discretized 
Brownian paths and along 40 individual paths with the SROCKC2-6 method with h = 
5x28. 


sider constructing methods with higher strong global convergence orders and 
better stability properties. 
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